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Abstract 

ESPRIT  is  a  successful  algorithm  for  determining  the  constant  directions  of  arrival  of  a  set 
of  narrowband  signals  on  an  array  of  sensors.  Unfortunately,  its  computational  burden  makes 
it  unsuitable  for  real  time  processing  of  signals  with  time- varying  directions  of  arrival.  In  this 
work  we  develop  a  new  implementation  of  ESPRIT  that  has  potential  for  real  time  processing.  It 
is  based  on  a  rank-revealing  URV  decomposition,  rather  than  the  eigendecomposition  or  singular 
value  decomposition  used  in  previous  ESPRIT  algorithms.  We  demonstrate  its  performance  on 
simulated  data  representing  both  constant  and  time-varying  signals.  We  find  that  the  URV- 
based  ESPRIT  algorithm  (total  least  squares  variant)  is  effective  for  time- varying  directions-of- 
arrival  using  either  rectangular  or  exponential  windowing  techniques  to  diminish  the  effects  of 
old  information. 
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1  Introduction 


The  esprit  algorithm  [8]  is  a  very  successful  means  of  determining  directions-of-arrival  (doa) 
of  a  set  of  narrowband  signals  impinging  on  an  array  of  sensors  with  translational  invariance.  It 
handles  array  geometries  almost  as  general  as  those  of  the  music  algorithm  [9]  at  a  significant 
computational  savings. 

A  key  limitation  of  both  the  MUSIC  and  ESPRIT  algorithms  is  that  they  require  0(nm2)  work 
to  process  n  data  samples  from  an  array  of  m  sensors.  Unfortunately,  at  the  heart  of  each 
algorithm  is  the  solution  of  an  eigenvalue  problem  or  singular  value  problem,  and  there  is  no  good 
way  to  economically  update  the  previous  information  in  order  to  reduce  this  computational  cost 
when  new  signal  information  arrives,  so  processing  one  new  data  sample  also  requires  0{nm?) 
or  0(m3)  work.  This  burden  makes  both  algorithms  unsuitable  for  real-time  computation.  Even 
using  parallel  processing  0(m2)  processors  are  required  to  achive  0{m)  time. 

Some  attempts  have  been  made  to  reduce  the  updating  complexity  by  maintaining  an  approx¬ 
imate  singular  value  decomposition  (e.g.,  [5]),  but  we  believe  that  better  results  can  be  obtained 
using  an  alternate  decomposition. 

Recently,  Stewart  [10]  has  introduced  anew  matrix  decomposition,  the  rank-revealing  URV  de¬ 
composition,  that  requires  somewhat  less  computation  than  either  a  singular  value  decomposition 
or  an  eigendecomposition  but  reveals  much  of  the  same  information.  A  key  advantage  of  the 
URV  decomposition  is  that  it  can  be  updated  in  time  0(m )  on  an  array  of  m  processors.  This 
means  that  algorithms  that  previously  depended  on  eigendecomposition  or  singular  value  de¬ 
composition  might  now  be  practical  for  real  time  applications  if  the  URV  decomposition  can  be 
successfully  substituted.  Boman,  Griffen,  and  Stewart  [1]  have  already  exploited  this  fact  in  ac¬ 
celerating  the  music  algorithm.  The  purpose  of  the  present  work  is  to  investigate  the  use  of  the 
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urv  algorithm  in  time-varying  signal  processing  using  esprit.  Roy  and  Kaliath  [8]  noted  that 
the  use  of  the  eigendecomposition  or  singular  value  decomposition  was  in  some  sense  unnecessary 
for  ESPRIT,  because  all  that  is  required  is  a  basis  for  the  range  space  of  a  certain  matrix.  They 
lacked  a  more  economical  matrix  decomposition  that  reliably  provided  this  information.  We  show 
in  this  work  that  the  rank- revealing  URV  decomposition  is  the  appropriate  tool. 

The  following  two  sections  give  brief  descriptions  of  the  total  least  squares  ESPRIT  algorithm 
and  the  URV  decomposition.  Section  4  discusses  the  algorithm  tv-esprit  ,  a  time- varying;  imple¬ 
mentation  of  the  ESPRIT  algorithm.  Experimental  results  are  summarized  in  §5,  and  computa¬ 
tional  details  of  the  urv  decomposition  are  given  in  an  appendix. 

2  The  ESPRIT  Algorithm 

Consider  d  narrow-band  plane  waves  simultaneously  incident  on  a  planar  array  of  to  sensors, 
arranged  in  m/2  doublet  pairs,  where  m  is  an  even  integer.  The  displacement  A  between  sensors 
in  a  pair  is  a  constant,  but  the  location  of  each  pair  is  arbitrary.  The  wave  sources  are  assumed 
to  be  located  in  the  same  plane,  and  the  location  of  each  source  is  specified  by  a  single  parameter 
6{  6  [0,  27t] ,  the  doa  of  the  ith  source. 

For  notational  convenience,  we  will  denote  data  related  to  the  first  sensor  in  each  pair  by  a 
subscript  X ,  and  data  related  to  the  second  by  Y ,  and  all  vectors  will  be  column  vectors. 

Given  data  from  the  array  of  sensors,  the  DOA  estimation  problem  is  to  locate  the  directions 
of  the  sources  of  radiating  energy  that  is  being  detected  by  the  sensors.  If  the  signals  are  assumed 
to  be  narrow-band  processes,  with  the  same  known  center  frequency  o>0,  then  the  DOA  problem 
can  be  described  by  a  simple  model:  the  relationship  between  the  unknown  signal  s(t)  £  Cd  and 
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the  sensor  output  zx(t)  6  £m/2  and  zY(t)  €  Cm /2  is  given  by 


*x(t) 

=  As(t)  +  cx(t), 

(1) 

zY{t) 

=  A$s(t)  +  eY(t), 

(2) 

z(t)  =  «(*)  +  e(t)  (3) 

\A*I 

where  e(f)  is  the  measurement  noise,  and  A  6  Cm/2*d  is  the  unknown  matrix  of  array  responses 
or  array  steering  vectors,  dependent  on  the  directions  of  arrival.  The  diagonal  matrix  $  is  also 
unknown,  and  is  related  to  the  phase  delays  between  the  sensors  in  each  doublet  pair: 


&•  =  ej"°Asln0l/c,  i  = 


Our  task  is  to  estimate  the  number  of  signals  d  and  the  directions-of-arrival,  9{,i  =  1,  •  •  *,d.  For 
this,  it  is  sufficient  to  estimate  the  matrix  4>,  and  this  is  the  underlying  idea  in  ESPRIT  . 


2.1  esprit  for  constant  source  directions 

The  esprit  algorithm  exploits  the  array  geometry  in  the  following  way.  In  the  absence  of  noise 
(i.e.,  i  =  Q),  the  range  of  the  matrix  A  is  the  same  as  the  range  of  the  matrix  Zx  formed  from 
n  columns  of  sensor  outputs  (n  sufficiently  large)  and  that  is  the  same  as  the  range  of  ZY ,  taken 
over  the  same  time  interval.  From  the  information  in  the  basis  vectors  for  the  range,  it  is  possible 
to  construct  a  matrix  that  is  similar  to  the  unknown  diagonal  matrix  $,  and  thus  determine  $ 
from  the  eigenvalues.  Determining  the  basis  vectors  and  computing  the  matrix  similar  to  $  can 
be  done  in  many  different  ways,  but  each  variant  of  ESPRIT  has  these  two  phases: 

A.  Obtain  a  basis  for  the  range  space  of  A. 
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B.  Compute  $  by  finding  the  eigenvalues  of  a  matrix  similar  to  #,  and  compute  the 
DOAs  from  this. 

Phase  A  contains  the  bulk  of  the  computational  work.  Two  different  approaches  have  been 
taken  to  determining  the  range  space  basis:  computing  a  singular  value  decomposition  of  the  data 
matrix  ZH  or  computing  an  eigendecomposition  of  the  estimated  covariance  matrix 

Rzz=-it^UH(t)  =  -zzH  ■  (5) 

For  comparison,  we  describe  one  previous  version  of  the  ESPRIT  algorithm  [8],  one  based  on 
the  singular  value  decomposition.  The  svd  of  the  matrix  ZH  is 

ZH  =  U  EFH, 

where  U  is  an  n  X  m  matrix  and  V  is  a  m  X  m  matrix,  both  having  orthonormal  columns.  The 
matrix  S  is  diagonal  containing  the  singular  values  in  descending  order:  oq  >  •  •  •  >  crm.  The  first 
d  columns  of  V  form  the  required  basis.  Note  that  U ,  which  has  the  same  dimensions  as  ZH ,  is 
not  needed  in  the  analysis  and  can  be  discarded. 

The  first  four  steps  correspond  to  Phase  A;  the  remaining  ones  process  the  range-space  basis. 
Our  URV-based  algorithm,  presented  in  §4,  differs  in  the  steps  marked  with  asterisks. 

Algorithm:  esprit(svd) 

1)  Obtain  the  data  measurements  Z. 

2*)  Compute  the  singular  value  decomposition  of  ZH : 

ZH  =  UY,Eh. 

3*)  Estimate  the  number  of  sources  d  (the  rank  of  ZH)  using  An  Information  Criterion 
(AIC)  or  the  Minimum  Description  Length  (MDL)  criterion. 
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4)  The  basis  for  the  signal  subspace  (the  range  of  ZH)  is  Ez,  equal  to  the  first  d 


columns  of  E. 

5*)  Partition  Ez  into  m/2  x  d  blocks  as 

Ex 
Ey 

and  compute  the  singular  value  decomposition  of  ( Ex,Ey ): 

(. Ex,Ey )  =  WVVH. 

6)  Partition  V  into  d  X  d  blocks  as 

Vu  V12 
V21  V22 

and  calculate  fa,  the  eigenvalues  of  — V12V22  • 

7)  Estimate  the  doas  from  fa  using  (4). 

Computing  the  SVD  of  the  data  matrix  becomes  impractical  if  the  matrix  is  too  large.  An 
alternate  approach  is  to  work  with  the  sample  covariance  matrix  Rzz,  which  is  m  X  m.  The 
eigenvectors  corresponding  to  the  d  largest  eigenvalues  of  Rzz  form  a  basis  for  the  approximate 
rangespace. 

A  different  approach  which  might  be  used  is  to  compute  a  rank-revealing  QR  factorization  of 

the  data  matrix  ZH  [3].  The  goal  of  this  algorithm  is  to  factor  ZH  as 

\ 

R\\  R\2 

0  R22  1 

where  R22  is  (m  —  d)  X  (m  —  d)  and  small  in  norm.  This  is  done  using  a  two-phase  procedure,  in 
which  the  standard  QR  algorithm  (perhaps  using  pivoting)  is  used  to  compute  an  initial  factor¬ 
ization.  Then  the  presence  of  a  small  singular  value  in  Rn  is  detected  by  condition  estimation 
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algorithms.  If  such  a  value  is  found,  the  dimension  of  R22  is  increased  and  a  permutation  and  up¬ 
date  of  R  is  performed  to  create  small  elements  and  move  them  to  i722-  The  condition  estimation 
process  is  repeated  as  necessary.  The  columns  of  Q  form  the  basis  for  the  range  space  of  ZH ,  but 
the  orthogonal  basis  for  the  range  of  Z  is  not  directly  available. 

Note  that  the  basis  for  the  null  space  of  a  matrix  is  not  really  needed  in  Step  6  of  the  algorithm. 
Since  VHV  —  I,  we  know  that  V^Vw  +  Vff V21  =  0,  and  —V^V^1  =  so  a  basis  for  the 

range  space  suffices. 

The  svd  ,  eigendecomposition,  and  QR  approaches  can  all  yield  reliable  estimates  of  the  range 
space,  but  they  have  drawbacks.  None  is  suitable  for  time-varying  computation,  since  none  of 
the  decompositions  can  be  updated  efficiently  as  new  data  arrive.  Parallelization  of  the  methods 
seems  to  require  an  array  of  0(m 2)  processors  in  order  to  achieve  0{m)  or  O(n)  time.  Recently, 
less  computationally  expensive  approximations  to  the  singular  value  decomposition  have  been 
proposed,  but  it  is  not  clear  how  they  behave  with  transient  data  [2,  5,  6,  7].  Moreover,  these 
approximations  still  require  0(m2)  processors  to  update  in  parallel. 

The  main  goal  of  our  work  is  to  develop  a  technique  for  Phase  A  that  is  similar  to  the  singular 
value  decomposition  in  reliability  but  faster,  more  parallel,  and  easier  to  update. 

2.2  Requirements  for  a  Time- Varying  esprit  algorithm 

In  the  time  dependent  problem,  the  sensors  receive  a  new  data  sample  at  each  time  unit.  There 
are  two  common  approaches  to  discounting  the  old  data  in  order  to  develop  reliable  estimates  of 
the  current  DOAs. 

The  first  is  the  (rectangular)  windowing  method.  In  this  method,  only  the  most  recent  n  data 
samples  are  retained,  and  earlier  ones  are  discarded  as  being  irrelevant.  One  way  to  express  this 
is  to  multiply  the  data  Zi,  the  data  collected  at  time  i,  by  a  rectangular  window  function  of  fixed 
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dimension  n  defined  as 


Wn,N(i )  =  < 


1  i  =  N,  N  —  —  n+1 


0  otherwise 

This  function  works  like  a  window  frame  that  only  admits  n  pieces  of  data,  and  we  shift  it  forward 
every  trial  to  use  the  most  recent  n  samples  of  data.  At  time  N,  we  work  with  the  data  samples 


(-2A-TI+1  5  —5  ZN-\,  Zn). 

The  second  approach  uses  forgetting  factors  to  discount  the  data  in  a  more  gradual  way.  As 
each  new  data  sample  is  received,  the  old  data  is  multiplied  by  a  number  //  between  0  and  1.  In 
this  method,  all  preceeding  data  is  multiplied  by  an  exponential  window  function,  so  at  time  N 
we  work  with  the  data  samples 

Since  the  dimension  of  this  matrix  is  growing,  it  is  essential  that  our  numerical  technique  works 
with  a  compressed  form  of  the  data,  a  matrix  whose  dimension  is  independent  of  N. 

Updates  to  the  the  range  space  cannot  be  computed  in  real  time  using  either  the  eigendecompo- 
sition  or  the  singular  value  decomposition,  so  we  turn  to  a  new  decomposition,  the  rank-revealing 
URV. 


3  The  urv  Decomposition 

The  determination  of  an  approximate  range  or  nullspace  for  an  n  x  rn  matrix  ZH  is  a  recurring 
problem  in  many  areas  of  signal  processing.  We  will  assume  that  ZH  is  a  matrix  consisting  of  true 
values  plus  measurement  errors,  and  that  in  the  absence  of  error,  ZH  would  have  rank  d  <  m. 
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The  rank  revealing  URV  decomposition  expresses  ZH  in  the  form 

(  \  (  u  \ 

R  F  v¥ 

ZH  =  ( Uz  Ux)  =  urvh, 

1°  a)VH 

where  the  columns  of  U  and  V  are  orthonormal,  R  and  G  are  upper  triangular  of  orders  d  and 
m  —  d,  and 

«  +  \\G\\%  <  e2. 


This  decomposition  reveals  that  ZH  is  within  e  of  a  matrix  of  rank  d  (distance  measured  in  the 
Frobenius  norm).  From  this  it  is  seen  that  ||Z^Vx||f  <  so  that  Vz  is  a  basis  for  the  approximate 
rangespace  of  the  true  (errorless)  data  matrix  and  V±  is  a  basis  for  the  approximate  null  space. 
Moreover,  it  is  not  necessary  to  carry  U  along  to  compute  and  update  the  decomposition,  so  that 
the  storage  requirements  are  modest. 

The  rank- revealing  urv  decomposition  is  not  unique;  the  singular  value  decomposition  and 
the  rank-revealing  QR  factorization  are  both  special  cases,  and  tools  such  as  plane  rotations  and 
condition  estimators  are  common  to  the  entire  family.  The  usefulness  of  the  URV  approach  hinges 
on  choices  that  make  it  rank- revealing  and  inexpensive  to  update,  requiring  only  the  use  of  a 
fixed  number  of  plane  rotations.  The  presence  of  the  matrix  V  (distinguishing  urv  from  qr) 
ensures  that  an  explicit  basis  for  the  approximate  nullspace  is  available.  The  flexibility  of  using  a 
triangular  factor  rather  than  a  diagonal  one  (distinguishing  urv  from  svd)  greatly  improves  the 
efficiency  of  the  update. 

The  urv  decomposition  has  the  following  properties  [10]: 


•  It  requires  0(m 2)  storage. 

•  It  furnishes  the  matrix  V±  explicitly. 

•  It  can  be  updated  in  0(m2)  time. 
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•  With  m  processors,  it  can  be  updated  in  0(m)  time. 


Given  a  rank  revealing  URV  procedure,  we  can  implement  an  ESPRIT  algorithm  based  either 
on  manipulation  of  ZH  or  of  Rzz-  Updating  the  factorization  when  new  data  arrive  proceeds 
in  three  steps,  described  in  detail  in  [10].  We  add  the  new  row  to  the  R  matrix  and  determine 
whether  the  matrices  F  or  G  have  become  too  large  in  norm  (based  on  a  tolerance  7).  If  so, 
we  tentatively  increase  our  estimate  of  the  rank  and  perform  some  rotations  that  change  R  and 
V.  We  then  reduce  R  to  upper  triangular  form  without  changing  V .  We  need  to  decide  whether 
R  has  become  rank  deficient;  this  is  the  standard  condition  estimation  problem.  A  tolerance  of 
7-7  is  used  to  test  for  possible  deficiency  here.  If  there  is  a  potential  rank  deficiency,  then  a 
pivoting  procedure  is  performed,  modifying  R  and  V,  to  make  the  last  column  of  R  small.  Then 
a  refinement  step  is  performed  in  order  to  bring  the  triangular  matrix  closer  to  diagonal  form.  It 
suffices  to  perform  the  rank-deficiency  check  at  most  twice  per  update. 

For  the  forgetting  factor  method,  the  current  matrix  ZH  is  replaced  by 


l  u  ^ 
pZH 


\ 


) 


where  p  is  the  forgetting  factor  and  z  is  the  new  data  vector.  Once  we  have  a  rank-revealing 
URV  decomposition  of  the  data  matrix  ZH,  it  is  an  0{m 2)  process  to  compute  the  rank-revealing 
URV  decomposition  of  the  modified  matrix,  We  expect  the  columns  of  F  and  G  to  have  size 
approximately 

I  (m  —  d ) 

7SV 

(measured  in  the  2-norm),  where  a w  is  the  standard  deviation  of  the  noise. 

If  the  windowing  method  is  used,  one  additional  step  is  needed.  We  first  downdate  the 
decomposition  by  deleting  the  first  data  row  in  ZH ,  and  then  perform  the  updating  algorithm 
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with  fi  =  1.  The  details  of  downdating  are  discussed  in  the  appendix.  Again,  the  total  work  is 
0(m2),  but  the  most  recent  n  data  samples  must  be  saved,  increasing  the  storage  requirements 
to  0(nm).  We  expect  the  columns  of  F  and  G  to  have  size  approximately 

7  =  \J(m  —  d)noN, 

so  again  one  should  choose  a  good  estimate  for  <7jv  for  the  rank  test  tolerance. 

4  The  tv-esprit  Algorithm 

Algorithm:  Time  Varying  esprit(urv)  Suppose  that  we  already  have  a  rank  revealing 
urv  decomposition  of  the  data  matrix  from  the  previous  time.  For  windowing,  we  also  save  the 
most  recent  n  data  samples. 

1)  Obtain  the  new  data  sample  z. 

2)  Update  the  previous  rank  revealing  URV  decomposition  by  downdating  and  updating 
the  previous  factors  if  windowing  is  used,  or  updating  the  previous  factors  if  forgetting 
factors  are  used. 

3)  Estimate  the  number  of  sources  d  (i.e.,  the  rank  of  R )  using  a  tolerance  of  r  times 
the  standard  deviation  of  the  noise.  (The  parameter  t  is  chosen  by  the  user.) 

4)  The  basis  for  the  signal  subspace  (the  range  of  ZH )  is  Ez,  equal  to  the  first  d 
columns  of  the  V  factor  in  the  URV  decomposition. 

5)  Partition  Ez  into  m/2  X  d  blocks  corresponding  to  the  X  and  Y  sensors: 

/ 

Ex 

Ez  = 

V Er 
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and  compute  a  rank- revealing  URV  decomposition  of  (Ex,  Ey): 

( Ex,Er )  =  W*VH. 

6)  Partition  V  into  d  x  d  blocks  as 

Vu  V12 

V21  V22 

and  calculate  <f>i,  the  eigenvalues  of  -VuV^1  (or,  equivalently,  the  eigenvalues  of 

VnHV£). 

7)  Estimate  the  doas  from  using  (4). 

In  Step  3,  we  chose  r  to  be  3  or  6.  In  Step  5,  we  set  r  =  1,  7  =  .5  if  the  signal  sources  were 
separated  by  at  least  10°,  and  r  =  10,  7  =  .05  to  force  more  refinement  of  the  subspace  in  case 
the  signals  were  close. 

The  URV  decomposition  in  Step  5  requires  0(d3)  time  since  we  will  get  a  different  Exy  every 
trial  and  can  not  apply  the  updating  technique  to  reduce  the  required  operations.  Thus  the  total 
time  per  step  is  0(m2  +  d3),  and  in  practice  the  number  of  signal  sources  d  is  often  much  less  than 
the  number  of  sensors  m.  Even  if  d  is  large,  Steps  1-3  require  only  0(m2),  and  thus  it  should  be 
possible  to  keep  up  with  the  incoming  data,  and  if  necessary  update  the  DOA  estimates  at  less 
frequent  intervals. 

In  comparison  with  the  SVD  algorithm,  the  basis  for  the  signal  subspace  changes  less  frequently: 
only  if  the  algorithm  believes  that  the  subspace  is  changing  dimension.  Thus  if  the  signals  are 
moving  slowly  compared  with  the  sampling  rate,  we  expect  the  URV-based  algorithm  to  require 
many  fewer  updates  of  the  doas  than  the  svd  algorithm. 
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5  Experimental  Results 


In  this  section,  we  will  present  some  simulation  results  that  illustrate  the  performance  of  the  two 
algorithms:  esprit(svd)  and  esprit(urv). 

We  use  a  five-pair  (m  =  10)  linear  array  with  pair  spacing  A/4.  The  pairs  are  placed  on 
a  line  of  length  4A  with  relative  location  [0, 1,3, 5.5, 7]  (A/2).  The  two  signals  are  narrow-band 
with  signal  to  noise  ratio  (SNR)  23dB  and  20dB,  respectively.  The  noise  is  Gaussian,  cind  the 
algorithms  were  tested  with  duplicate  data  samples  in  order  to  make  a  fair  comparison. 

We  say  that  an  algorithm  failed  at  a  particular  time  if  it  estimated  more  or  fewer  than  two 
signals. 

The  first  example  concerned  two  fixed  signal  sources  located  at  24°  and  29°  and  with  50% 
initial  correlation.  For  each  trial,  we  estimated  the  DOAs  based  on  100  data  samples,  and  we  ran 
2000  trials.  Figure  6.1  shows  a  histogram  and  tabular  summary  of  the  results.  Both  algorithms 
were  quite  successful.  This  shows  that  we  are  not  sacrificing  much  accuracy  by  substituting  the 
more  economical  urv  for  the  svn  . 

The  second  example  tracked  two  signals  with  time- varying  doa  using  an  exponential  forgetting 
factor  of  0.9.  We  considered  two  cases: 

•  close  sources:  periodic  doa  of  10°  +  5°  sin(27rn/360)  and  20°  +  5°  sin(27rn/240), 

•  well-separated  sources:  -10°  +  10°  sin(27rn/360)  and  40°  +  5°  sin(27rn/240), 

where  n  =  1,3,..  .,719.  We  ran  100  trials  for  each  case.  The  average  error  was  defined  to  be  the 
mean  of  the  differences  between  estimated  and  actual  DOAs. 

Figure  2  and  3  show  that  the  DOA  estimations  from  the  two  algorithms  are  quite  similar.  The 
error  histogram  in  Figures  4  and  5  show  that  the  variance  of  the  relative  error  between  actual 
and  estimated  DOAs  for  esprit(urv)  is  quite  small.  When  n  £  [375,425],  the  DOA  s  of  the  close 
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sources  are  almost  equal,  and  both  algorithms  were  fooled  into  believing  that  there  was  a  single 
source.  In  the  figures,  we  plotted  the  results  from  the  most  recent  successful  trial,  so  we  have 
straight  lines  for  both  tracks  in  this  interval.  Note  that  the  svd  based  algorithm  is  not  a  practical 
computational  tool,  since  the  data  matrix  grows  far  too  large,  but  it  provides  a  useful  standard 
of  comparison  to  assess  the  quality  of  the  URV-based  algorithm.  Figure  6  shows  the  number  of 
floating  point  operations  as  a  function  of  the  number  of  samples  in  a  single  trial.  The  work  for 
the  URV-based  algorithm  grows  as  m2,  while  the  work  for  the  SVD-based  algorithm  grows  as  bib2. 

As  a  third  example,  we  used  the  same  data  as  in  the  second,  but  a  rectangular  window  of  size 
20  rather  than  forgetting.  The  window  shifted  2  points  for  each  trial.  Again,  the  results  given 
in  Figure  7  and  8  and  error  histograms  given  in  Figure  9  and  10  show  that  the  estimations  are 
acceptable.  Figure  11  plots  the  number  of  floating  point  operations  as  a  function  of  the  number 
of  samples.  Again  it  is  clear  that  the  work  for  the  URV-based  algorithm  is  quadratic,  while  that 
for  the  SVD  algorithm  grows  cubicly. 

As  a  final  example,  to  demonstrate  tracking  of  instantaneously  changing  signals,  we  assumed 
that  there  were  two  signal  sources  located  at  24°  and  29°,  each  with  SNR  23  dB,  and  that  the 
signals  alternatively  appear  and  disappear.  We  took  a  rectangular  window  of  size  10.  Figure  12 
shows  the  similar  performance  of  the  two  algorithms. 

These  experimental  results  lead  us  to  believe  that  the  URV-based  ESPRIT  algorithm  can  be 
successfully  used  for  real-time  tracking  of  time- varying  signals. 

6  Summary 

We  have  presented  a  new  variant  of  the  esprit  algorithm  that  has  potential  for  real-time  tracking 

of  moving  signals.  It  lias  the  following  features: 
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•  The  storage  requirement  is  0(m2)  (plus  mn  for  rectangular  windowing). 

•  The  work  per  update  is  0{m 2  +  d3). 

•  It  performs  nearly  as  well  as  svD-based  algorithms,  at  greatly  reduces  computational  cost, 
and  admits  an  efficient  parallel  realization. 
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6.1  Appendix:  Downdating  the  URV 


The  computations  related  to  updating  the  URV  decomposition  when  new  data  appear  have  been 
described  in  [10],  so  we  concentrate  in  this  section  on  the  downdating  algorithm  when  data  are 
deleted.  This  downdating  algorithm  takes  0(m2)  time  and  computes  the  urv  factorization  when 
the  first  data  row  in  the  window  is  deleted. 

Suppose  that  we  already  have  a  URV  factorization  of  an  n  X  m  data  matrix  ZH  with  n  >  m. 
Specifically,  only  the  R  and  V  matrices  are  known.  Then 

(r\ 


ZH  =  U 


v0/ 


v, 


where  U  is  an  unknown  n  x  n  unitary  matrix.  Then 


/ 


ZHVH  = 


z»VH 


\  (  H\  (  -  \ 
Hi 


R 


(6) 


\ZHVH  j 


/ 
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where  zx  and  ux  are  the  first  rows  of  ZH  and  U  respectively  and  ZH  has  rank  d.  Now,  we  want 
a  URV  factorization  of  the  matrix  Z^ .  Referring  to  (  6),  our  problem  is  equivalent  to  downdating 
a  QR  factorization  of  ZHVH  [4]. 

First,  we  compute  the  vector  .  The  first  m  components  of  this  vector  are  determined  by 
solving  the  linear  system 

w?R  =  gVH.  (7) 

Although  the  block  G  of  R  is  small,  the  corresponding  subvector  of  zxVH  is  also  small,  and  the 
process  is  well  determined. 

Columns  m  +  1  through  n  of  U  can  be  chosen  to  be  any  basis  for  the  null  space  of  ZH ,  and 
without  loss  of  generality,  we  assume  that  columns  m  +  2  through  m  have  a  zero  in  the  first  entry. 
Therefore,  the  (m  -f  l)st  entry  of  ux  is  p,  =  1  -  || || ,  and  the  remaining  ones  are  zero. 

Now,  we  determine  a  sequence  of  left  rotations  that  rotate  ux  into  a  multiple  of 

the  first  unit  vector.  Rotation  Jk  modifies  the  kth  and  ( k  +  l)st  rows  in  order  to  zero  out  the 
( k  +  l)st  element.  That  is 


Ji--'Jmu1  =  aex  (8) 

Applying  these  rotations  to  the  triangular  factor  in  (6)  from  the  left  will  yield  a  Hessenberg  matrix 
H:  a  matrix  with  zeros  below  the  first  subdiagonal.  We  now  have  the  relation 


/ 


\  /  T  \ 
a  0T 


V 


J 


y  o  u 


H  . 


(9) 


z”VH 
ZH  VH 

Note  that  the  first  column  of  the  unitary  matrix  must  be  a  multiple  of  the  first  unit  vector  since 
|a|  =  1,  which  forces  all  of  the  other  entries  in  the  column  to  be  zero.  Therefore  we  have  a 
decomposition  of  the  matrix  ZH  VH  as  a  matrix  with  orthonormal  columns  times  a  triangular 
matrix  R  consisting  of  rows  2  through  m  +  1  of  H,  and  our  task  is  completed. 
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It  requires  about  4m  multiplications  and  2m  additions  to  perform  one  rotation.  Since  we 
apply  m  rotations  in  the  downdating  processing,  the  total  time  required  is  0(m2). 

In  order  to  preserve  the  rank- revealing  properties  of  the  URV  decomposition,  we  now  need  to 
check  whether  the  triangular  matrix  has  dropped  in  rank,  a  procedure  described  in  [10]. 
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Figure  1:  Histogram, results  means  and  variances  of  ESPRIT(svd)  and  ESPRIT(URV)  for  fixed 
sources 
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Figure  2:  Estimated  time-varying  DOAs  for  close  sources  using  exponential  windowing. 
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Figure  3:  Error  histogram  for  DOAs  for  close  sources  using  exponential  windowing. 
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Figure  4:  Estimated  time- varying  doas  for  well- separated  sources  using  exponential  windowing. 
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Figure  7:  Estimated  time- varying  doas  for  dose  sources  using  rectangular  windowing. 


17 


moving  DOA  (degree)  moving  DOA  (degree) 


50 

40 

30 

20 

10 

0 

-10 

-20 

-30 

0  100  200  300  400  500  600  700 

time 

50 

40 

30 

20 

10 

0 

-10 

-20 

-30 

0  100  200  300  400  500  600  700 

time 

Figure  9:  Estimated  time- varying  doas  for  well-separated  sources  using  rectangular  windowing. 
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